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BOUNDARY AND INTERFACE CONDITIONS FOR HIGH ORDER FINITE 
DIFFERENCE METHODS APPLIED TO THE EULER AND NAVIER-STOKES 

EQUATIONS 

JAN NORDSTROM * AND MARK H. CARPENTER t 

Abstract. Boundary and interface conditions for high order finite difference methods applied to the 
constant coefficient Euler and Navier-Stokes equations are derived. The boundary conditions lead to strict 
and strong stability. The interface conditions are stable and conservative even if the finite difference operators 
and mesh sizes vary from domain to domain. Numerical experiments show that the new conditions also lead 
to good results for the corresponding nonlinear problems. 

Key words, high-order finite- difference, numerical stability, interface conditions, summation-by-parts 

Subject classification. Applied and Numerical Mathematics 

1. Introduction. In many computational problems, low order finite difference methods (second order 
or less) are not accurate enough. Examples in which high-frequency components in the solution must be 
resolved by using high order finite difference methods (HOFDM) include aeroacoustics, turbulence and 
transition simulations, the propagation and scattering of electromagnetic waves, and simulation of reactive 
flows at high speeds [1], [2], [3], [4], [5], [6]. The efficiency [7] of HOFDM can be used either to increase 
the accuracy for a fixed number of mesh points or to reduce the computational cost for a given accuracy by 
reducing the number of mesh points. 

The main reason that low order finite difference methods are used in practical calculations is because of 
the difficulty that arises for HOFDM near the boundaries of the computational domain. On a Cartesian mesh, 
it is quite easy to derive nonsymmetric boundary operators that have high formal accuracy; the difficulty is 
to derive highly accurate and stable operators. In [8] and [9] HOFDM are constructed based on the work in 
[10] and [11]. In these strictly stable schemes, the growth rates of the analytic and semidiscrete solution are 
identical. Strict stability is obtained by constructing discrete operators that satisfy a summation-by-parts 
(SBP) rule which mimics the integration- by- parts rule in the continuous case. For calculations over long 
times, strict stability is very important beacuse it prevents error growth in time for fixed Ax. 

In [12] it was shown that many G-K-S stable [13] (convergence to true true solution as Ax — >• 0) scalar 
schemes were not strictly stable. Moreover, many scalar schemes that were both G-K-S stable and strictly 
stable exhibit time growth when they are applied to systems of equations. The underlying reason for the 
error growth in time caused by the way the mathematical boundary conditions were imposed. An orthogonal 
projection operator is used to impose the mathematical boundary conditions in [14] and [15]. In the so called 
SAT (simultaneous approximation term) procedure [16], a linear combination of the boundary conditions in 
the form of a forcing function and the differential equations is solved near the boundary. Both these methods 
impose the correct boundary conditions and preserve strict stability. 
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Administration under NASA Contract No. NAS1-97046 while this author was visiting the Institute for Computer Applications 
in Science and Engineering (ICASE), Mail Stop 403, NASA Langley Research Center, Hampton, VA 23681-2199. 
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Another important concept is strong stability. An approximation is strongly stable if the solution, 
including the values at the boundary points, can be estimated in terms of all data in the problem [17]. The 
stability estimate in the strongly stable case leads directly to the error estimate if no extra or numerical 
boundary conditions are necessary. Stability analysis using the Laplace transform technique leads to strong 
stability if the Kreiss condition is satisfied; see [18] and [9] paper IV. Note that strict stability leads to strong 
stability, but strong stability does not imply strict stability. 

Most investigations regarding HOFDM are done on linear hyperbolic model equations with constant 
coefficients on a uniform mesh. However, nonlinear Navier-Stokes calculations on nonuniform meshes have 
been performed [19]. One of the conclusions in [19] was that the treatment of the metric derivatives is a 
crucial point for nonsmooth meshes. This problem is analyzed in [20] where so called mimetic dif ference 
operators (discrete operators with the same symmetry properties as the continuous operators) are derived. 
In [15], strict stability for parabolic and hyperbolic systems in curvilinear coordinates on a single domain 
were investigated. 

Generating a smooth grid around a complex configuration can be very difficult, if not impossible, and is 
often the most time-consuming aspect of the solution procedure. This fact has limited the use of HOFDM in 
practical calculations to the small class of simple geometries which can be smoothly mapped onto the unit 
cube. In this paper we consider a structured multiblock approach in which each subdomain is discretized by 
using a discrete operator with the SBP property. The subdomains are patched together to a global domain 
by using suitable interface conditions. This technique was used in [21], [22] and [23] for Chebyshev spectral 
methods. 

In [24] , stable and conservative interface conditions for HOFDM applied to the scalar advection-diffusion 
equation on multiple domains were derived. In each subdomain the step size was constant but significantly 
different from that in the adjacent subdomains. Also, the finite difference operators could vary from sub- 
domain to subdomain. In this paper we will generalize the results in [24] and extend the analysis to the 
one- dimensional constant coefficient Euler and Navier-Stokes equations. 

The rest of this paper will proceed as follows. In section 2, some basic definitions are given. In section 3, 
the Navier-Stokes equations on conservative, primitive, and characteristic variable form are given. In section 
4, the continuous problem is analyzed, while the discrete problem is investigated in section 5. Numerical 
experiments are performed in Section 6 and we summarize and draw conclusions in section 7. 

2. Definitions. Consider the linear initial boundary value problem 

wt — Pw-\-5F{x,t) ,x£Ct ,t> 0, 

(2-1) w = Sf(x ) ,xGd , £ = 0, 

Lew — Sg(t ) , x e T , t > 0, 

where P is the differential operator and Lc is the boundary operator. The initial function <5/, the forcing 
function SF, and the boundary data Sg are the data of the problem; w denotes the difference between a 
solution with data f,F,g and one with data f -\-5f,F+5f,g + 6g. There are many concepts of well posedness, 
see [17]. Here we consider the following definition. 

Definition 1. The problem (2.1) is strongly well posed if the solution w is unique, exists, and satisfies 

( 2 - 2 ) IMIq + IMIr* < K^iWSfWl + fqSFWf, + IMIr)*}, 

Jo Jo 

where K c and r) c may not depend on SF, Sf,Sg. || • |[n and || • [|r are suitable continuous norms. 
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The semidiscrete version of (2.1) is 

(■ u>j)t = Qwj+5Fj(t) ,Xj G 0, ,t> 0, 

(2.3) Wj — 5fj ,Xj ef) ,t = 0, , 

LdWj = Sg(t) ,Xj e T , £>0, 

where Q is the difference operator approximating the differential operator P, 5Fj is the forcing function, Sfj 
the initial function, Lo the discrete boundary operator where numerical boundary conditions are included, 
and Sg the boundary data. It is assumed that (2.3) is a consistent approximation of (2.1). 

Closely related to the concept of well posedness is the concept of stability. 

Definition 2. The problem (2.3) is strongly stable, if for a sufficiently fine mesh, the solution Wj 
satisfies 

(2-4) IMft + f IMIr<« < K^iWSfWl + f(\\SF\\l + HMIrW, 

Jo Jo 

where Kd and rjd may not depend on SFj, Sfj,Sg. j| • j|^ and || • ||r are suitable discrete norms. 

Definition 3. The approximation (2.3) of (2.1) is strictly stable if the analytical and discrete growth 
rates (see (2.2) and (2.4)) satisfy 


(2.5) 


Vd < r]c + 0{ Ax), 


where Ax is the mesh size. 

For later reference we also define some useful matrix operations; see [25]. 

Definition 4. Let A be a p x q matrix, B be an m xn matrix, and Ii the l x l identity matrix, then 


ao,oB 

ao,q-iB ^ 

, Ii®B = 

/BO • 
0 B 

0 \ 
0 

ip-i,o£ • 

Q-p— l,q— 1£ / 


\ 0 0 

.. B J 


The pxq block matrix A®B and the l xl block diagonal matrix Ii<S>B are called Kronecker products. There 
are a number of rules for Kronecker products (see [25]). In this paper we will make use of, 


( 2 . 6 ) 


{A ® B)(C & D) — {AC) <g> ( BD ), {A ® B) T = A T 0 B T . 


The following lemma will be used frequently below; it is a direct consequence of the first rule in (2.6). 

Lemma 1 . Let A be an mxm matrix, B be annxn matrix, A = I n ®A and B — B<S>I m , then AB = BA. 
Proof : The first condition in (2.6) leads to AB — {I n ® A){B <g> J m ) = B ® A = {B ® / m )(/ n 0 ^4) = BA. 
□ 


3. The Euler and Navier-Stokes equations. The one-dimensional constant coefficient Navier- 
Stokes equations in primitive (W), characteristic (C), and conservative (Q) variable form are 

(3.1) W t + AW X = eBW xx , C t + AC X = eXC xx , Q t + F(( = eF^ 

respectively. With e = 0, equation (3.1) becomes the one-dimensional constant coefficient Euler equations. 
The overbar is used to denote variables with a constant state. The relation between W, C, Q where W = 
(p, u, T) t is 

(3.2) C = RSW, Q = TW 
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where 



( 

— l/v^T 

l/y/2 

— v/7-1/ 2 7 ’ 

R, = 

V7- 1/7 

0 



V 

1/V27 

1/V2 

^7-1/27 / 


1 

( d 2 / 1/7 

0 

0 ) 

S = y/2 

0 

pc 

0 


I 

l 0 

0 p! \ j 7(7 — / 


/ 1 

u 

\ c 2 /(t(7- 1)) +w 2 /2 


0 0 

P 0 

pu p/( 7(7 - 1 )M^) 


Note that i?i? T = 73 . 

The transformation (3.2) implies that the matrices and fluxes in (3.1) are 


(3.3) 


(3.4) 


(3.5) 


(3.6) 


(3.7) 


A = 


( u 

c 2 /lP 

V 0 


p 

u 

( 7 -l )c 2 M^ 


° ^ 

1/7 Ml , 

u ) 


( 00 0 

0 (A + 2p)/p 0 

0 0 7 k/{Prp) 


A = (RS)A(RS)- 1 = 

X = (RS)B(RS)- 1 = i 
F 1 = TAW = TAT “ 1 Q, 


u — c 

0 

0 ^ 

0 

u 

0 

0 

u 

+ c ) 


a(f) 

6-4> 

CK 0 

a 2 4> 

—a4> 


—a<p 

6 + j> 


P v = T5W, = TBT~ X Q X . 


The dependent variables and parameters p,u,T,p,c, Moo,^,A,/c,Pr ,7 and e are respectively the density, 
x,y,z components of the velocity, the temperature, the pressure, the speed of sound, the free-stream Mach 
number, the shear and second viscosity, the coefficient of heat conduction, the Prandtl number, the ratio of 
specific heats, and the inverse Reynolds number. The notations 6 — (A + 2/2) /p, 4> = (7 — 1 )R/(Prp),ct — 
— 1 ) has also been introduced. 

4. The continuous problem. In this paper we will consider interface conditions between subdomains. 
However, interface conditions are closely related to boundary conditions; therefore, we start with the single 
domain problem. 
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4.1. The continuous single domain problem. To make the presentation self-contained, some results 
in [27] are included in this section. Consider the Navier-Stokes equations on characteristic form, 


C t + AC x = eXC xx + F(x,t) 

,t> 0 

, — 1 < x < 1, 

C = f{x) 

,t = 0 

,-l < x < 1, 

L-\C = g~i(t) 

,t > 0 

, X — 1, 

L + iC — g+i(t) 

,t > 0 

,X = +1, 


where C — ( pcu — p, a(pc 2 —p), pcu + p) T ,0 < e << 1 and L- 1 , L+\ are the boundary operators. For u > 0, 
there is inflow at x — — 1 and outflow at x = 1. 


4.1.1. Well posedness. Let 

<**-/>* (U,U) = \\Ut, ||t/||? = \U\l = _, + |£/g =+ i 

denote the L 2 scalar product, the L 2 norm, and the boundary norm respectively. The energy method applied 
to (4.1) leads to 

lie'll? = [ C t AC - 2 eC T XC x }lz+l - 2e(C x ,XC x ) + 2(C,F). 

The boundary conditions (see [27] and [22]) 

(4-2) L^C= (±+M C - eXC x = g„ u 


(4-3) L+iC ~ | (A 2 |A|) C-^C i | = {g+i}i, *=1,2, 

where |Aj = diag{ ]A X |, |A 2 |, |A 3 |) leads to 


||C||, 2 = -2e(C il ,XC I ) + 2(C,F) 

(4-4) - [C t A,C - 2C r S _ll«.-i - [C T \ 0 C + 2C T g +1 ] x=+u 

where g +1 = ( 91 , 52,91 - (2 /a)g 2 ) T and 


(4.5) 


Aj = |A|, A 0 = 


! | Ai| 0 (|Ai| — Ai)/2 

0 |A 2 | 0 

V (I Ai. | — Ai)/2 0 | A 3 1 


Integration of (4.4) leads to 

pT T 

lie'll 2 + e’ T { 2 e / ( C x ,XC x )e~'> t dt + | [ ||C|| 2 e-”*dt} 

Jo 1 Jo 

( 4 - 6 ) < e ” T {||/|| 2 + ? [ T Mh^dt + - [ T ||F|| 2 e-"'di}, 

0 Jo V Jo 

where 0 < p < 1, 6 — min \dj\, D = JA|17, and 


H = diag(H u 1,H 3 ), 


TT = 1^3 | ~ |Al[ 

1^3 | -F | A 3 [ ’ 


#3 


|A 3 1 - |Ai| 
I A 3 1 + |Ai| 


Note that (4. 2), (4. 3) reduce to the characteristic boundary conditions for the Euler equations as e — ■» 0. 

Uniqueness follows directly from the estimate (4.6). Existence can be shown by using the Laplace- 
transform technique or via difference approximations; see [26] and [28]. Since (4.6) is of the form (2.2), we 
can conclude that the following theorem holds. 

Theorem 1. The problem (4-1) with the boundary conditions (4-2), (4-3) is strongly well posed. 
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4.2. The continuous multiple domain problem. In this section we split the domain [—1,1] into 
[—1,0] and [0, 1] and focus on the interface problem at x = 0. The two coupled problems are 


U t + AU X 

— eXU xx + F(x , t) 

,t> 0 

,-l < x < 0. 

U 

= f( x ) 

,t = 0 

,-l < x < 0, 

L-iU 

= 9-i(t) 

,t > 0 

,x = -1, 

L 0 (U-V) 

= 0 

,t> 0 

,x = 0, 

V t +AV x 

= eXV xx -f- F(x , t ) 

,t> 0 

, 0 < x < +1, 

V 

= /(*) 

,t = 0 

,0 < x < +1, 

Lo(V-U) 

= 0 

o 

A! 

-w 

,x = 0, 

L+iV 

= 9+i(t) 

,t> 0 

+ 

II 


respectively. The characteristic variables in the left [—1,0] and right [0, +1] domain are U and V respectively. 
The coupling between (4.7) and (4.8) is given by the operator Lq. 

By subtracting (4.1) from (4.7-4.8), by transforming the problem on [0, +1] onto [—1,0] via the trans- 
formation x — > — £, and finally by replacing £ with x, we obtain 


(4.9) 


where 


and 

(4.10) 


+ A i[) x 

c Xip xx 

,t> 0 

,-l < x < 0 


= 0 

II 

o 

,-l < x < 0, 

L-Jj 

= 0 

,t> 0 

i—l 

1 

II 

L+iV 

= 0 

,t > 0 

,x = -1, 

Lo(U-V) 

= 0 

,t > 0 

II 

o 


i '> = 


u 

V 


u -c 
v-c 


A = 


-f A 0 
0 -A 


X = 


X 0 
0 X 


LiU = 


(A + |A|) ~ 


U-eXU x ,L +1 V = 


(A-|A|) 


V + eXK 


A i 

X f ? 1 

) i 


1 , 2 . 


4.2.1. Well posedness. The energy method applied to (4.9) leads to 


Mi = - 2e'ip T Xip x ]* =0 1 - 2e(-0 x , Xip x ). 


The analysis of the single domain problem implies that the boundary terms at x = — 1 are negative semidef- 
inite with the boundary operators (4.10). At the interface x = 0, we have 

V> T A - 2eijj T X'ip x ] x=z o = 



( ) 


/ 0 A -eX 0 N 


( XJ-V > 

1 

u + v 


A 0 0 -eX 


U + V 

2 

(U - V) x 


-eX 0 0 0 


(U-V) x 


{ (U + V)* ) 


^ 0 -eX 0 0 j 


\ ( u + v) x / 


Well posedness for the Euler equations (e = 0) requires U — V = 0 since A is nonsingular. With that 
choice we get 

[^ T Mj - 2e^F X iJ) x \ x=z q = -2eU T X(U + V) x = -2 e{{RS) T U) T B(W 1 + W 2 ) x , 
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where ( RS)~ 1 U = Wr, ( RS)~ 1 V — Wr denotes the primitive variables in the left and right domain respec- 
tively. The structure of B (see (3.4)) and a transformation to the original coordinate system lead to the 
following theorem. 

Theorem 2. If Theorem 1 holds and the interface conditions 


(4.12) 



U-V 

(U-V) x 


— 0, D i 


1 0 1 
1 a -1 


are used, then (4-7) and (4-8) are strongly well posed. 

Remark. The problems (4.7) and (4.8) are strongly well posed in the sense that the solutions can be 
estimated in ternms of the data in the corresponding one domain problem (4.1). 

Remark. The condition (4.12) in primitive variable formulation is 


h 

eD 2 


W l -W r 
(W l - W R ) X 


= 0, D 2 


0 1 0 \ 

0 0 1 J ' 


5. The discrete problem. Let UfDU be the numerical approximations of the scalar quantities u and 
u x respectively. The approximation VU of the first derivative 


VU = P~ 1 QU, Pu x — Qu — PT el , |T e i| = Q{ Ax m ,Ax n ) 


satisfies the SBP rule 

(5.1) (U,VV) P = U N V N -U 0 V 0 -(VU,V)p 
where 

(5.2) ( U,V)p=U T PV ; P=P T , Q + Q t = D, D = diag[- 1,0, ..,0,1] 

and 0 < p m i n AxI < P < p max AxI. Operators of the SBP type arise naturally with centered difference 
approximations; for examples see [11], [29], [12], [30]. 

The second derivative can be obtained by applying the first derivative operator twice. Such an ap- 
proximation satisfies the SBP rule (5.1) exactly. However, there are drawbacks with such a procedure. A 
second derivative formed in that way is unnecessarily wide and inaccurate and can lead to odd-even mode 
decoupling. A second derivative operator with the following properties, 

(5.3) V 2 U = P~ 1 RU, Pu xx — Ru — PT e2 , T e2 = <D{ Ax m ,Ax n ), 

(5.4) R = (S t M + D)S, 

was suggested in [24]. The matrix D is given in (5.2); M is positive definite, i.e., U T MU > 0 and 0 < 
m min AxI <M< m max AxI. 

S is a diagonal matrix with a discrete representation of the first derivative on the first and last rows, 
{Su} 0 = {Vu} o = u x (x 0 , t ) + T e3 , {5u} n = { T>u} n = u x (x n , t) + T e3 , 
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where |T e3 | = 0(Ax r ) and 



The second derivative defined in (5.3) and (5.4) satisfies a modified SBP rule We have 

{U,V 2 V) P = U n {VV} n - U 0 {VV } o - ( SU) t M{SV ). 

The notation |T e i|,|T e 2 | = 0(Ax m , Ax n ) and |T e 3 | = 0(Ax r ) means that the approximation of the 
differential operator is accurate to order m in the interior of the domain, to order n at the boundary and 
that the approximation of the boundary conditions is accurate to order r. The relation between the different 
orders of accuracy, i.e., m,n,r is discussed in section 5.1.2 below. 

Examples of first and second derrivative approximations are given in (A.1)~(A.5) in appendix A. The 
approximations are second order accurate in the interiour of the domain and first order accurate at the 
boundary. This means that for (A.1)-(A.5) we have m = 2 and n = 1. 

So far we have considered difference approximations of scalar quantities. The corresponding approxima- 
tions for vector quantities are defined by using Kronecker products (see definition 4). The spatial operators 
V, V 2 and the matrices that define them are of the form B 0 J 3 in this paper. As an example, P~ l Q means 
(P -1 (g) h){Q 0 I 3 } — P~ l Q 0 h- In the sequel, that notation is implied. 

Let H — H T > 0; for later reference we introduce the notations 

(5.5) (U,V) h = U t HV, (U,U) H = \\Uf H , rilLH^?=o + |E/iL, 

5.1. The discrete single domain problem. We introduce a uniform mesh Xi = — 1 +iAx,xo — 
— l,x n = +1. The finite difference approximation of (4.1) with the SAT technique [16] for boundary condi- 
tions is 

C t +AVC = eXV 2 C + F 

(5.6) + P~ 1 {a- 1 (L£ 1 C - 0_i)e_i +< t +1 (L£ 1 C - g +1 )e +1 ], 

C( 0) = /, 

where 

(5.7) V = P-'Q ® I 3 , V 2 = P~ 1 R®I 3 , 

(5.8) R - QP-'Q <g> h or R = (S T M + D)S 0 / 3 , 

(5.9) A = / n 0) A, X = I n ®X, e_i = (l,...0 ) T 0/ 3 , e +1 = (0, ...1) T 0 J 3 . 

The unknown diagonal matrices er_ 1 and cr + 1 will be determined below. 
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5.1.1. Stability. The energy method leads to 

^j|C||! = -C T (AQ + Q T A)C + eC T (XR + R T X)C + 2(C,F) P 
(5-10) + 2C% - g_i] + 2C%a +1 {L° 1 C - g +1 ]. 

The definition of the first derivative operator P~ 1 Q and Lemma 1 leads to 
(5-11) -C T {AQ + Q T A)C = C^ACo - C%AC n . 

The definition of the second derivative operators R — (— S T M -fl D)S and R = QP~ l Q yields 

C t (XR + R t X)C = -2C 0 XVC 0 + 2 C n XVC n 

(5.12) - {SC) t {XM + (XM) t ){SC) 

C t (XR + R t X)C = -2C 0 XVC 0 + 2C n XVC n 

(5.13) - 2 (P- 1 QC) T PXP~ 1 QC, 

respectively. 

By introducing (5. 11), (5. 12) and (5.13) into (5.10) we get 

d 

J t \\ c fp + 2 c{DC, XVC) H = [C T KC - 2 eC T XVC)tz° n + 2 (C,F) P 

+ 2Cg , <7_i[L£ 1 C- 9 _ 1 ] 

(5-14) + 2Cfla + i[L± 1 C — g+i], 

where the scalar products and norms are defined in (5.5) and 

R = QP~ 1 Q => H = P, 

R = (- S t M + D)S^H = (S(P- 1 Q)- 1 ) t ( M + 2 MT XSiP^Q)- 1 ). 

The boundary operators L I ? 1 ,L+ 1 are the discrete versions of (4.2)-(4.3), with one important modification. 
In [27] it is shown that the two outflow conditions in (4.3) determine the value of the last row of XC x in 
terms of the in-going characteristic variable and boundary data; i.e., (4.3) implies that 

(5-15) {-eXC x } 3 = - Al ~ |Al| C! + gl - (2 /a)g 2 , x = +1. 

To exphcitly incorporate (5.15) into (5.6) we use 

(5.16) £giC={ (A+ 2 |A|) C — eXVcj = g. u 

(5-17) ^°C=[ (A ~ |A|) C- e XPe) =g +1 , 

where (< 7 + 1)3 is equal to the right-handside of (5.15). The boundary conditions (5. 16), (5. 17) inserted in 

(5.14) yields 

~\\Cfp = -2 e(VC,XVC) H + 2(C,F)p 

+ {C T [+A + a i(A + |A|)]C} i=0 + {C r [-2f.X - 2err_ t X\VC},,„ 

+ {C T [-A + <7 +1 (A - |A|)]C} i= „ + {C T [+ 2eX - 2ca n X\VC}, =n 

(518) + {4, C'C^Xi - |Ai|)} i= „ + 2C 0 Vi - 2 Cj s+1 . 
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The choice. 


(5.19) (j-i = -I 3 , (7 + i = / 3 , 

leads to 


(5.20) 


l|C||? = - 2 e(VC,XVC)„ + 2 (C, F)p 

- [C t A/C - 2C T 5_i]i =0 - [C t A 0 C + 2CVi]i=», 


i.e., a growth rate which is exactly the same as in the continuous case (compare (5.20) with (4.4)). The 
definitions of A/, Ao are given in (4.5). Integration of (5.20) leads to 

||Cj|f, + e r?£>T {2e f ( VC,XVC) H e-™ t dt+ ^ f \\C\\l dt) 

Jo 2 Jo 

(5.21) <e I?DT {||/||p + f [ T \\g\\l D e^dt + 4- f \\Ff P e-^di). 

od Jo Vd Jo 

The estimate (5.21) is similar to (2.4) and hence (5.6) is a strongly stable approximation. The problem 
(5.6) is also strictly stable (we can choose t)d — V and Sd = see (4.6), (2. 5)). We can summarize the result 
in the following way. 

Theorem 3. The approximation (5.6) of the problem (4-1) is both strictly and strongly stable if (5.19) 
holds. 


5.1.2. Accuracy. The problem describing the deviation Ej = C(xj,t) — Cj(t) between the exact 
continuous solution and the discrete approximation given by (5.6) is 


E t + A VE = eXV 2 E + T 

(5-22) + p-^*_ 1 (L» 1 E)e-i+cr + i(L» 1 E)e+ 1 }, 

£7(0) = 0. 

T — T do + T bc is the truncation error. T DO and T BC comes from the approximation of the differential 
operator and the approximation of the boundary conditions respectively. The truncation errors have the 
general structure 


(5.23) 


rpDO 


/ 0(Ax n ) \ 
0{ Ax m ) 

0( Ax m ) 
v 0( Ax n ) ) 


T bc = 


( 0( Ax^- 1 )) \ 

0 


{ 0( Ax( r ~V J 


In [31] and [32] it is shown that difference approximations to mixed hyperbolic-parabolic equations 
retain the accuracy of the interior scheme (0(Ax m )) if a finite number of points (independent of the total 
number) are closed with boundary stencils (0(Ax n )) that are one order less accurate. A requirement for 
that conclusion is that an energy estimate holds, which in turn means that the mathematical boundary 
conditions must be approximated to the order of the internal scheme. The discussion above implies that 
that n = m — 1 and r — m is necessary. 

We will now apply the theory in [31] and [32] to the type of difference approximations considered in this 
paper, i.e., where difference operators of the SBP type are used together with a penalty formulation for the 
boundary conditions. 
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First, we split E and the T into two parts, i.e., E = E 1 + E 2 and T = T 1 + T 2 where 

f o \ ( 90 \ 


(5.24) 


T 1 


9 1 

9n— 1 

V o J 


0{ Ax m ), T 2 = 


0 

V 9n J 


= 0( Ax^^). 


Next, we use the energy method to estimate E 1 . The energy method applied to (5.22) with E, T replaced 
by £J 1 ,T 1 , and the conditions (5.19) leads directly to 

H^Hp < G( Ax m ). 

Finally we use the Laplace-transform technique to take care of the boundary error and estimate E 2 . 
So far, the treatment has been general. However, in order to keep the algebraic complexity at a reasonable 
level, we now need to simplify and be specific. We will consider the inviscid (e = 0) Euler equations at an 
inflow boundary approximated with the second order scheme given by (A.l) and (A.2) in appendix A. The 
half-plane problem obtained by Laplace-transforming (5.22) with E,T replaced by E 2 , T 2 becomes 

sE 2 +A(E 2 -E$) = cr_i(A + |A|)l?o + Axgo, 

( 5 - 25 ) «je? + A(JS? +1 -.E?_ 1 )/2 = 0, j> 1, 

0, j 


E i 


oo 


where s — sAx. The second and third equations in (5.25) lead to 

/i\ , /OX o \ 

0 + <T2 1 "b a 3 0 

0 / V 0 / 1 / 


(5.26) 


E 2 = a, 


j 

v 3’ 


_ (®/ Aj) + \/l + ( s/Xj ) 2 , Aj > 0 

(5.27) Kj = 0 , Aj = 0 , 

— v / i + (W , Aj < o 

where the branch of the square root is the one with positive real part for Re(s ) > 0. The case when A j = 0 
presents no problem; it only reduces the number of equations in (5.6). In the sequel, we assume A j 0. 
The first equation in (5.25) leads to 

(5.28) E{s)a = Ar vg 0 , E(s ) - diag(s + XjKj - (Aj + crLfyA j + j Aj j)), j = 1, 2, 3. 

A nonsingular E(s)), i.e., 

(5.29) det(E(s)) ± 0, Re(s) > 0, 
and (5. 26), (5. 27) lead to 

\E 2 \ < const . | Axgo |, Re(s ) >0, j > 0; 

i.e., the Kreiss condition is satisfied. Parsevabs relation and the fact that E 2 (t) cannot depend on go(T) for 
t <T leads to 

r-t /»t 


f \E 2 \ 2 dt < const, f |Ax^o | 2 ^5 j > 0 
Jo Jo 
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Fig. 5.1. The mesh close to the interface at x = 0. 


and finally, since go — O(Ax), 


||E 2 ||p < 0( Ax 2 ). 

It still must be shown that (5.29) holds. The inviscid condition for strict stability Xj +<ri 1 (Aj + 1 Ay | ) < 0 
(see (5.18) and (5.27)) which implies s+XjKj = |Aj|-y/l + ( s/Xj ) 2 > 0, leads directly to (5.29). The procedure 
to estimate the boundary error at an outflow boundary is exactly the same as in the inflow case. We can 
summarize the result in the following Theorem. 

Theorem 4. The approximation (5.6) of the problem (4-1) with e = 0 is second order accurate if 
Theorem 3 holds and the first derivative operator T> = P~ X Q is given by (A.l) and (A. 2) in appendix A. 

Remark. The procedure that was exemplified above to prove accuracy in the second order accurate case 
is general. The last step where one uses the Laplace-transform technique to estimate the boundary error 
E 2 is not necessary i) if the boundary stencils have the same order of accuracy as the internal stencil, i.e., 
n — m and ii) if the approximation of the mathematical boundary conditions is one order more accurate, 
i.e., r = m + 1. 


5.2. The discrete multiple domain problem. A finite difference approximation of the coupled 
problems (4.7) and (4.8) is 

eXV 2 L U + F + BT 0 

~ Vo) + *Z((Z > L U) n - ( V R V)o))e L 

f 

eXV 2 R V + F + BT m 

Pr^r{V o ~ U n ) + * v R {{V R V)o - (: V L U) n ))e R 

/• 

The characteristic variables in the left (subscript L) [xo — — l,x n — 0] and right (subscript R) [x 0 = 0,x m = 
+1] domains are U and V respectively, see Figure 5.1. BTo,BT m denote the boundary terms at x = ±1 
respectively. Definitions of X>,£> 2 , A, X,e_i,e + i are given in (5. 7), (5. 8), (5. 9), and ez, = (0, , , 1) T <g> / 3 , e R = 
(1„,0 

The values of er_i and a + i that lead to strict and strong stability for the discrete single domain 
problem are given in (5.19). We must still determine , a R , cr R . Note that the difference operators 

can be different in the left and right domains and that Axl / A xr. 


(5.30) 


U t + kV L U = 
+ 

C/(0) = 
Vt+AVnV = 
+ 

V(0) = 


5.2.1. Conservation. To calculate the strength and speed of a shock with finite mesh size, one needs 
a conservative scheme. Let us start by considering a continuous problem in conservation form, ut + f x — 
0,|a;| < 1, £ > 0. Integration over the domain leads to 

d f +1 

— y ^ u dx + f+i /—i — 0, 
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i.e., the total change of u in the domain is only due to the flux through the boundaries. Note that integration 
of f x over the the domain reverses the differentiation process and leaves information only at the boundaries. 
Let F,T>F denote the numerical approximations of /, f x . The discrete SBP derivative satisfies 

(5.31) f x — Vf = T e i , £>/ — P~ 1 Qf, T e i = G(Ax r ). 

Multiplying (5.31) with the operator l T P where l T = [1, 1, ..., , 1] <g> I p (/ has p components) and observing 
that f + i — f- 1 = f x dx leads to 

r+l 

l T Pfx = J fx dx + 0(Ax r ). 

The operator l T P is the discrete integration operator. This operator reverses the process of differentiation, 
leaves information only at the boundaries, and converges to the continuous integration operator as Ax — > 0. 
We can now prove the following theorem. 

Theorem 5. The approximation (5.30) of the problem (4-1) is conservative if 

(5.32) 4 — A = 0, o v L -o\ + & = 0, 

where the matrices A and X are given in (3. 5), (3. 6). 

Proof : Multiplying (5.30) with i fPr and IrPr leads to 

(l lPlU + llPnV), = -{IIQl'KU + llQ R AV) + t(l T L R L XU + I t r RrXV) 

+(*L - °rWn - Vo) + (?Z - al)(VU N - Wo) 

(5.33) +2 (U, F) Pl + 2(V, F)p h + 

where BT includes the boundary terms at x — ±1. To obtain (5.33) we have made use of Lemma 1. 

The inviscid terms can be written 

(5.34) IIQl'AU + llQ R k = -(AU)o + (kU) n - (AV) 0 + (AV) m . 

Next, we consider the viscous terms. Both R = QP~ 1 Q and R = (— S T M + D)S lead to 

I t l QlPPQlXU + iIQrPr'QrXV = - {XVU ) o + (XVV)n 

(5.35) - (XVV)o + (XVV) m . 

By inserting (5.34) and (5.35) into (5.33), neglecting the boundary terms at x = ±1, letting F = 0, and ap- 
plying condition (5.32), we obtain (IJ^PlU + l])PnV)t = 0; i.e., the approximation (5.30) is conservative. □ 

5.2.2. Stability. We start with the following observation. 

Remark. Stability of the one domain problem does not imply stability of the multiple domain problem. 
Stability means that the solution can be estimated in terms of the (bounded) boundary data. In a multiple 
domain problem, the boundary data are made up of the solution(s) in the other domain(s). Boundedness of 
the data would require an a priori assumption. 

The main result of this paper is given below. 

Theorem 6. The approximation (5.30) of the problem (4-1) is both strictly and strongly stable if 

(5.36) <rI=<reX, o‘ L = \(k- f3eX - 5h), + 6 > 0, 

2 ZOLR ZOiL 
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and if Theorem 3 and 5 hold. oil,&r denotes the minimal eigenvalue of P if R — QP~^Q and the minimal 
eigenvalue of (M + M T )( 2 if R= (~~S T M + D)S. The matrices A,X are given in (3. 5), (3.6). 

Proof : Strict and strong stability of (5.30) follows if the interface treatment at x = 0 is of a dissipative 
nature. For that reason we neglect the terms at the boundaries x — ±1 and use F = 0. The energy method 
leads to 

j t (\Mk + lir||p R ) = -u t (aq l + q t l K)u - v t (aq r + q t r a)v 

+cU t (XR l + R t l X)U + eV T (XIt n + R‘„X)V 
+2UZ(al(U n - Vo) + c r\{VU n - VV 0 )) 

(5-37) +2V 0 t (4(V 0 - E/„) + al(VV 0 - VU„)). 

Equations (5. 11), (5. 12) and (5.13) lead to 

(5-38) -U t (AQl + QlA)U = UfAUo - U%AU n 

(5.39) ~V t (AQ r + Q t r A)V = VfAVo - V%AV m 

(5-40) U T (XR L + R t l X)U < -2U 0 XVU a + 2U„XVU n - 2a L VUl XVU n 

(5.41) V t (XR r + R t r X)V < -2V„XVVo + 2V m XVV m - 2a R VV? XVV 0 . 

By inserting (5.38)-(5.41) into (5.37) and neglecting boundary terms at x = ±1 we obtain 

gait'll l, + imilj < w t ew, 

where 



U n 

\ 


( 2<t£ — A 

-(?l + °r) 

*Y + eX 

^>4 

b 

1 

\ 


Vo 


,E = 

-(*l + °i) 

2 a R + A 


*Y ~eX 



vu n 


*Y + eX 


—2 OiL,eX 

0 



Wo 

) 


b 

1 

*%-eX 

0 

—2 aReX 

/ 


The problem (5.30) is strictly and strongly stable if E is negative semidefinite. E is an almost full 
matrix; to obtain explicit stability conditions, simplifications of E are necessary. The energy method applied 
to the continuous multiple domain problem leads to (4.11) which suggests that the variables 



( U n - Vo 

\ 


+1 

-I 

0 

0 > 

1 

U n + Vo 


s-4_ 

+1 

+/ 

0 

0 

vt 

VU n - Wo 


’ V2 

0 

0 

+1 

-I 


\ VUn + Wo 

J 


0 

0 

-f -I 

+ 1 J 


are of interest. The use of these variables and the conservation conditions in Theorem 5 leads to 

2*1 + eX eX \ 

0 0 

— {otL + olr)cX (olr — a.L,)eX 
(cxr — (xl)cX — (o:l + Q-R)eX y 

To show that E\ is negative semidefinite, first assume that the first condition in (5.36) holds. Secondly, 
add and subtract the matrix — 2/3eX to the upper left block in E \ . The condition for negative semidefiniteness 
(see Lemma 1) becomes 

(5.42) y\ (2(2 *{ - A) + 2/3eX) yi + e[(Y <g> R) T y 2 ] T [A E2 ® B][{Y ® R) T y 2 ] < 0, 


Ei = SES t = 


( 2(2*1 -A) 0 
0 0 
2 *Y + eX 0 


V 


eX 


0 
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where B — R T XR, Ar 2 — Y T E 2 Y and 


E 2 — 


-2 P 

(1 + 2a) 
1 


(1 + 2a) 
~{OLL + OCR) 
OCR ~ OiL 


&R ~ &L 
— {OiL + &r) 


The first term in (5.42) is nonpositive if the second relation in (5.36) holds. Negative definiteness that implies 
A e 2 < 0 is obtained if the third relation in (5.36) holds. □ 


5.2.3. Accuracy. In this section we will consider the accuracy close to the interface. The procedure is 
similar to the one used in section 5.1.2 for the single domain problem. The problem describing the deviations 
Uj = U (xj , t) — Uj (t) and Vj = V (x 3 ,t) — Vj (t) between the exact continuous solutions and the discrete 
approximations given by (5.30) is 


(5.43) 


u t + kv L u 

= 

(XV\U + t l 


+ 

n - Vo)+aY((V L U) n - (V n V)o))e L 

(7(0) 

= 

0 

V, + kv R v 

= 

tXV 2 R V + T R 


+ 

Pr(°'r(V -U) + <({PrV) 0 - ( T>LU) n ))e R 

F(0) 

= 

0. 


For simplicity, we have used the notation U — U and V = V. Note also that the terms at the boundaries 
x = ±1 are neglected. The treatment at the boundaries x = ±1 has been discussed in section 5.1.2. 

Tl and Tr are the truncation errors from the approximation of the differential operator and the interface 
conditions. The truncation errors have the general structure 



( '■ ) 


( o( AzSr 1 ') 

\ 

T l = 

O(Axf) 

5 Tr — 

Oi Axl) 



{ OiAx^-V) ) 


\ '■ 

) 


The discussion in section 5.1.2 on the size of the truncation error is applicable also for the interface problem. 

Following the procedure in section 5.1.2, see [31], [32], one splits up the errors in two parts, the first 
part (T£,Tft) contains the truncation error of the internal scheme, and the second part contains a boundary 
contribution (T^,T^) with one order lower accuracy. The structure of these errors are 


1 

/ 


\ i 

( : 

\ 

II 

1? 

0{ Ax?) 


11 

CM 



I 

\ 0 

) 

1 1 

< 

( o( a4”*-‘>) 

) 



( 0 



( o(a4" _1) ) 

\ 

Tr — 

Oi Ax%) 


^ T 2 r = 

0 



\ ^ 

) 


l : 

J 


Also the error is divided into two parts; i.e, we consider U = U 1 + U 2 and V — V 1 + V 2 . 

By using the energy method, U 1 and V 1 will be bounded by and T^. This procedure is straightfor- 
ward, entirely similar to the one in section 5.1.2 and will therefore not be repeated here. Suffice it to say 
that the stability conditions given in Theorem 6 lead to 


Uplift + II^IK < o(Ax n L ) + o(A*E). 
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To bound U 2 and V 2 in terms of T 2 and T 2 t requires use of the Laplace-transform technique. That analysis 
is given in detail below. 

Also in this case, we keep the algebraic complexity down by considering the inviscid (e — 0) Euler 
equations approximated with the second order accurate approximation given by (A.l) and (A.2) in appendix 
A. The problem for U 2 — U and V 2 = V obtained by Laplace-transforming (5.43) becomes 


(5.44) 


sjJJn T A(U n U n ~ i) 
+ A (Vi — Vb) 

SL&j + A(Uj + i — Uj - 1 )/2 

*RVj + MV Hl -V^)/2 
U 3 

Vj 


2cr I L (U n - Vb) + A x L g L 
2(T^(V r 0 - U n ) + A x R g R 
0, j < n - 1 

0, j > 1, 

0, j -> -oo 
0, j -> oo 


where s L = sAx L , s R = sAx R , g L = G( Ax^ 1} ) and g R = 0( Ax^ X) ). 
The last four equations in (5.44) lead to 





(°\ 

/ 0 

(5.45) Uj = <x| 

0 1 

(4r n +4 

i (4r n +4 

0 


o ) 


l 0 / 

1 1 


(4) 


3 \ j-n 


(5.46) 


Vj = a 


( 1 
0 

V o 


(4) j + 4 


0 \ 

i ( n 2 R y + <j 3 r 


0 
0 

1 


(4 y, 


Ai) \/l + (s/Xj ) 2 , \j > 0 

(5.47) 4 = 0 , Xj = 0 , 

-(s/Xj) + y/1 + ( s/Xj ) 2 , Aj < 0 


~{s/\j) + y/l + (s/Xj) 2 , Xj > 0 

(5.48) 4=0 , Xj = 0 , 

-(«Ai) - 4 1 + ( s Aj) 2 Aj < o 

where the branch of the square root is the one with a positive real part for Re(s) > 0. Also, in the case in 
which Xj = 0 presents no problem, only the number of equations in (5.6) is reduced. We assume Xj ^ 0 in 
the following. 

The coefficients <Jl = (cr\, 4? 4) T an< ^ a R = (4> 4> 4) T be determined by the first two equations 
in (5.44). They, together with the first condition in Theorem 5, lead to 


E(§i, s R ) 


VL 

VR 


( A x L g L 
\ A x R g R 


(5.49) 


E(sl, s r ) 


s R I 3 - Aac^ 1 + A - 2<t£ 

2 (4 - A ) 


2a[ 


s L h + A k r + A - 2 al 


16 


A nonsingular E(sl, sr ) leads via the Kreiss condition and Parseval’s relation; (see section 5.1.2) to the 
estimate 


\\U 2 \\ Pl < 0( Axl) + 0(Ax%), II^Hp, < 0(Ax|) + 0(Ax%). 

It still must be shown that (5.29) for E, defined in (5.49), holds. A direct calculation using (5.49), (5. 47), 
and (5.48) leads to 

3 

Det(E) = I] Gj = | A j | 2 ( 1 + + fc/A^l + (S fi /A,)2) 

3 = 1 

+ + {s L /\j) 2 yJl + ( sr/Xj ) 2 . 

Let a /1 + (§l/ \ j) 2 ~ tjl + ifr and yj\ + {sr/ \ j) 2 — t)r + i£r where r)i,,r]R are non- negative. A simple 
algebraic test reveals that the imaginary part and the real part of Gj cannot be zero at the same time if 
the inviscid condition for stability A — 2a I L < 0 in Theorem 6 holds. We can summarize the result in the 
following Theorem. 

Theorem 7. The approximation (5.30) of the problem (4-1) with e = 0 is second order accurate; i.e., 

\\U\\ Pl + \\V\\ Pr <0(AxI) + 0(Ax 2 r ), 


if Theorem 6 holds and the first derivative operator T> — P~ l Q is given by (A.l) and (A. 2) in appendix A. 

Remark. Also in the interface case, see section 5.1.2, the procedure to prove accuracy, which was 
exemplified above in the second order accurate case, is general. The last step in which one uses the Laplace- 
transform technique to estimate the errors XJ 2 andV 2 is not necessary i) if the stencils adjacent to the interface 
have the same order of accuracy as the internal stencil, and ii ) if the interface conditions are one order more 
accurate. 


5.2.4. The discrete multiple domain problem in conservation form. The discrete multiple 
domain problem (5.30) can be transformed to conservative form by multiplying the equations with I n +i ® 
T(RS)~ 1 ,I m+1 <gi T(RS)- 1 , respectively. The result is 


(5.50) 


Ut + P^lQLF 1 -eR c F v ) = 
Vt + Pn^QnF 1 - ( R r F v ) = 


+ (1 /2)P P [(Fl-F^) 

+ (l + 2a)e(FZ -FX)-FE] 

- (1/2)VP£--Fi) 

- (1 + 2 v)e(FX-FX)+F%], 


where F T = F 1 — eF v and 


Ff = ( SI 3 + ef3TBT~ l )(U n - V 0 ). Ff = ( SI 3 + e0TBT^)(V o - U„). 


In (5.50), the forcing terms and the boundary conditions at x ± 1 are neglected. 

6. Numerical experiments. By making one- dimensional computations using the nonlinear Euler and 
Navier- Stokes equations, we can check whether the theoretical conclusions drawn from the analysis of the 
constant coefficient problem agree with the results obtained in practice. 

In the calculations below, we use the second order scheme (given in (A.1)-(A.5)) and the fourth and 
sixth order schemes reported in [24]. To integrate in time, a five-stage fourth-order RK scheme [33] has been 
used. Consider the stability condition (5.36). In the calculations below we have used cr = —1/2 and the 
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Fig. 6.1. L 2 Errors in calculations using the Euler equations. 

conservative estimate a J L — A/2 — SI where S is determined through tests. Often we use S — 1.0. Equation 
(5.32) has been used to determine the other parameters. 

First, we consider a sound propagation problem. The computational results, obtained using the nonlinear 
Euler equations at Mach number 0.5, are compared with an exact solution of the linearized problem. In 
Figure 6.1 The errors for second, fourth, and sixth order schemes using one domain (lDom), four uniform 
domains (4Dom), and eight randomly spaced domains (Rand) are shown. Clearly, the order of accuracy is 
independent of the presence and location of the interfaces. Due to the small amplitudes (oc 10“ 7 ) used in 
the sixth order cases, we encounter round off, which can be seen as the kink on the sixth order results. 

Next, we consider a viscous shock propagation problem at Mach number 2.0 and Reynolds number 150. 
The exact solution of the Navier- Stokes equation for this case can be found in [34]. In Figure 6.2, the errors 
for second, fourth, and sixth order schemes using eight uniform domains (Unif) and eight randomly spaced 
domains (NonU) are shown. Also in this case, the order of accuracy is independent of the location of the 
interfaces. 

The curves in the sixth order case are not straight, see Figure 6.2. The curves are formed as a mean 
value of 15 simulations where different wave speeds ws from —0.25 to 0.5 are used. The individual results 
for each wave speed are given in Tables 6.1 and 6.2. Note: ws = 0 is stationary shock, and we have subsonic 
wave speeds for ws < 0.3. The results from uniform grid calculations are shown in Table 6.1; results from 
nonuniform grid calculations are given in Table 6.2. The convergence rate between the two grids is listed. 
The asymptotic limit approaches —6. Note that the trends are identical between the nonuniform and uniform 
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(*Y)] ai-*»iaanj 

Viscous Shock Propagation 



Fig. 6.2. L 2 Errors in calculations using the Navier-Stokes equations. 


cases. 

In Figure 6.3, the propagating shock ( ws = 0.25) for four different times is shown. In this case, the sixth 
order scheme and 24 gridpoints were used in each domain. 

Finally we will discuss two additional questions concerning accuracy and stability/efficiency. To investi- 
gate the influence of interface conditions on accuracy, we made the calculations illustrated in Table 6.3. The 
calculations are run to a physical time T — 3 at Mach number 2.0 and Reynolds number Re = 250. The 
sixth order SBP scheme is used, and the number of total points is 289 evenly distributed on the interval 
— 1/2 < x < 1. The parameter in the study is the number of subdomains, keeping the total number of 
intervals constant. The number of subdomains ranges from 1 to 24. For the case of 24 subdomains, the 
spatial operator involves 12 boundary stencils (fifth-order) and one sixth-order interior stencil. No further 
divisions are possible when using the sixth-order SBP operator. Note that this case is only marginally less 
accurate than the single domain case, for which the most points are discretized with sixth-order stencils. 

The previous study indicates that there is little loss of accuracy when subdividing the domain. There 
are, however, other costs associated with domain subdivision. Introduction of additional interfaces into the 
domain changes the resulting eigenspectrum of the semidiscrete operator. In [22], a reduction in the effective 
CFL, when using a penalty boundary procedure, was observed. We experience a similar reduction in the 
stability envelop as the number of subdomains is increased. 

In Table 6.4, a study compares the effective CFL of a singledomain calculation, with those from a 
comparable grid divided into eight subdomains. Plotted are the errors, and the maximum stable CFL as a 
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wave 


speed 


a: 96 128 

b:128 192 


192 256 

256 384 


384 

512 


-0.2500 

- 0.2000 

-0.1500 

- 0.1000 

-0.0500 

0.0000 

0.0500 

0.1000 

0.1500 

0.2000 

0.2500 

0.3000 

0.3500 

0.4000 

0.4500 

0.5000 


-4.4460 

-3.1743 

-1.5344 

-3.4447 

-4.7040 

-3.3093 

-5.6203 

-2.0256 

-4.9470 

-7.5646 

-4.7062 

-5.9644 

-4.9922 

-3.8538 

-0.7633 

-4.8735 


-4.4858 

-4.3626 

-6.5479 

-5.0035 

-4.8410 

-6.0257 

-3.0687 

-6.6065 

-5.1708 

-5.5734 

-3.0715 

-6.8890 

-5.0159 

-5.9798 

-2.4389 

-5.9665 


-4.8856 

-4.6413 

-2.8759 

-6.5253 

-5.0488 

-5.0116 

-3.8456 

-6.2051 

-5.5378 

-6.0670 

-3.4963 

-6.1815 

-5.5773 

-5.3015 

-5.8286 

-5.7257 


-5.2375 

-4.7217 

-7.4052 

-5.6487 

-5.4075 

-5.4064 

-6.8111 

-5.6289 

-5.4018 

-5.6418 

-4.1048 

-6.1886 

-5.4589 

-6.3515 

-4.9988 

-5.8033 


-5.5610 

-5.2825 

-5.6489 

-6.2999 

-5.6641 

-5.6608 

-6.4191 

-5.2130 

-5.7503 

-5.8898 

-5.8041 

-6.2484 

-5.1898 

-5.8163 

-7.4170 

-5.8160 


Table 6.1 


UNIFORM grid, 8 subdomains, refinements between grids a:b, 6th order explicit; CFL = 0.2. 


wave 

speed 

a: 96 
b:128 

128 

192 

192 

256 

256 

384 

384 

512 

-0.2500 

-3.9508 

-5.0473 

-5.6119 

-5.9785 

-4.5714 

-0.2000 

-3.8816 

-4.6038 

-4.8316 

-6.3334 

-6.9001 

-0.1500 

-4.4373 

-4.8932 

-5.1011 

-5.4514 

-5.8566 

-0.1000 

-7.6431 

-3.4782 

-7.2908 

-5.8925 

-5.9310 

-0.0500 

-3.2678 

-7.4661 

-3.7127 

-6.8562 

-6.2791 

0.0000 

-4.8406 

-4.7810 

-5.2543 

-5.5640 

-5.7606 

0.0500 

-7.8667 

-2.9978 

-7.8219 

-4.5670 

-4.7133 

0.1000 

-4.2532 

0.6385 

-2.8840 

-3.9126 

-6.1144 

0.1500 

-1.4577 

-5.1589 

-6.7711 

-5.0271 

-5.2246 

0.2000 

-4.2245 

-4.6373 

-4.5406 

-5.2829 

-5.9054 

0.2500 

-2.9734 

-4.0155 

-5.4352 

-5.4569 

-5.6145 

0.3000 

-4.2383 

-3.2435 

-4.3861 

-5.5030 

-5.8943 

0.3500 

-4.1902 

-3.5835 

-3.6667 

-4.6596 

-5.5589 

0.4000 

-3.4505 

-3.4125 

-4.8703 

-5.3783 

-6.3374 

0.4500 

-2.7380 

-2.9597 

-3.6868 

-4.5673 

-4.8980 

0.5000 

-4.1279 

-3.7378 

-3.8465 

-4.8352 

-6.1756 


Table 6.2 

NONUNIFORM grid, 8 subdomains generated randomly, refinements between grids a:b, 6th ordex exp LtexL^. CFL 
max/min ratio is 6-47. 
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Fig. 6.3. Viscous shock propagation, a domain with randomly spaced interfaces. 


Subdomains 

LOGioerror 

1 

-4.527 

2 

-4.584 

4 

-4.457 

8 

-4.643 

12 

-4.313 

16 

-4.467 

18 

-4.342 

24 

-4.358 


Table 6.3 


Variation of L 2 error on number of subdomains with grid density constant. 


function of Reynolds number for the two cases. Note that while the errors are nearly equivalent for the two 
test cases, the maximum CFL for the single domain case is nearly a factor of two larger. 

7. Summary and conclusions. We have analyzed boundary conditions and interface conditions for 
the one- dimensional Euler and Navier-Stokes equations. Both the continuous and semi-discrete problems 
have been considered. 

We have considered summation-by-parts operators and derived strictly and strongly stable boundary 
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Re 

LOGioerror 

GF L max 

LOGioerror 

CFLmax 


-2.154 

0.55 

DNC 


900 

-2.242 

0.55 

-2.265 

0.30 

800 

-2.347 

0.55 

-2.376 

0.30 

700 

-2.477 

0.60 

-2.517 

0.30 

600 

-2.637 

0.60 

-2.698 

0.30 

500 

-2.841 

0.60 

-2.935 

0.30 

300 

-3.429 

0.65 

-3.617 

0.30 

200 

-4.027 

0.65 

-4.185 

0.35 

100 

-5.741 


-5.699 

0.35 

40 

-7.892 


-7.331 

0.20 

20 

-9.535 


-8.637 

0.20 

10 



-10.665 

0.18 


Table 6.4 

Variation of CFL number and L2 error with Reynolds number for single and multiple domain cases. 


and interface conditions for the Euler and Navier-Stokes equations. We have also considered the question 
of accuracy, both in the general case and more specifically for a second order accurate approximation of the 
Euler equations. 

The interface conditions are stable and conservative even if the finite difference operators and mesh sizes 
vary from domain to domain. Numerical experiments which include a sound propagating problem and a 
viscous shock propagating problem show that the new conditions lead to accurate and stable results for the 
corresponding nonlinear problems also. 

It was also shown by numerical experiments that there is little loss of accuracy associated with domain 
subdivision. However, the introduction of interfaces into the domain changed the eigenspectrum of the 
semidiscrete operator and caused a reduction of the CFL number by approximately a factor of two. 


Appendix A. Stencils. We now present a few examples of the specific form of the stencils that have 
the SBP property. For more accurate stencils, see [24]. The second order accurate discretization matrix that 
approximates the first derivative T> = P~ 1 Q is 


(A.l) 


'-2 2 
-10 1 



-1 0 1 
-2 2 
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where 


(A- 2) 


' 1 
2 


‘-I l 

1 


-10 1 




1 


-10 1 

1 

2 - 


— 1 i_ 


The second order accurate discretization matrix that approximates the second derivative T> 2 — P 1 (— S T M+ 


D)S is 


(A.3) 


1 -2 1 

1 -2 1 


V = 


(A x) 4 


L 


1 -2 1 
1-2 1 


where 


(A.4) 


and 


(A.5) 



The matrix M can be shown to be positive definite (and symmetric). 
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